In part based on Fortran code from Furio Ercolessi and the CHEM126 course of Kalju Khan.
In [1]:
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())
Out[1]:
In lab 1 we worked with the Lennard-Jones potential. Other simple potentials include the harmonic potential
$$ \begin{equation} V = \frac{1}{2} k \left( r - r_{\text{eq}} \right)^2 \end{equation} $$with harmonic constant $k$ and equilibrium position $r_{\text{eq}}$, the Kratzer potential
$$ \begin{equation} V = D_0 \left( \frac{r-r_{\text{eq}}}{r} \right)^2 \end{equation} $$with dissociation energy $D_0$, and the Morse potential
$$ \begin{equation} V = D_0 \left[ 1 - e^{-\alpha (r - r_{\text{eq}})} \right]^2 \end{equation} $$where $\alpha$ is the Morse parameter.
For example, carbon monoxide (CO) consists of one carbon atom (atomic weight 12.011) and one oxygen atom (atomic weight 15.999) separated by 1.1283 Angstroms. This can be modelled by
Using standard minimization algorithms (e.g. scipy.optimize.minimize) show that, for suitable initial guesses, all of the above potentials have minima at the expected location.
By plotting the potentials, understand why the success of the algorithms changes.
Briefly test how changing the potential changes the number of steps required by the algorithm.
Also compute the force $\frac{dV}{dr}$ at the minimum point for each potential to check that it vanishes as expected.
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
rcParams['figure.figsize'] = (12,6)
from scipy.optimize import minimize
In [2]:
Using your velocity-Verlet integrator from lab 1, evolve the CO molecule using the Morse potential. Note that, unlike lab 1, the two particles now have different mass. It is perhaps easiest to scale it so that the carbon particle has mass $1$ and the oxygen mass $1.332029$. Do not use periodic boundaries! Try using $\Delta t = 0.001$ up to $t=1$. Start from the equilibrium position and show that it does not evolve. The move the location of the oxygen atom slightly and see how it evolves.
In [2]:
The paper of Praprotnik, Janezic, and Mavri suggests a potential for a single water molecule (one oxygen atom and two hydrogen with atomic weights $1.008$) of
$$ \begin{equation} V = \sum_{l = 1}^2 D_0 \left[ 1 - e^{\alpha \Delta r_{OH_l}} \right]^2 + \frac{1}{2} k_{\theta} \Delta r_{HH}^2 + k_{r \theta} \Delta r_{HH} \left( \Delta r_{OH_1} + \Delta r_{OH_2} \right) + k_{rr} \Delta r_{OH_1} \Delta r_{OH_2}. \end{equation} $$The first term is the Morse potential again, and the other terms are the intra-molecule bond potentials. $\Delta r_{OH_l} = r_{OH_l} - r_{OH_{\text{eq}}}$ is the stretch in the distance between the oxygen atom and the $l^{\text{th}}$ hydrogen atom $r_{OH_l}$ and its equilibrium value, and $\Delta r_{HH} = r_{HH_l} - r_{HH_{\text{eq}}}$ is the stretch in the distance between the hydrogen atoms. The parameter values are
Fix the oxygen atom at the origin and start the hydrogen atoms at $(\pm 0.8, 0.6, 0)$. Use the minimization techniques to find the equilibrium position of the atoms, comparing against the known structure of water.
NOTE: if using scipy's minimization routine, you may have to increase the tolerance as high as $10^{-4}$ to make it converge.
In [2]:
Using your integrator, evolve the molecule and see how it behaves.
In [2]: